Implement number theoretic transform for large integer multiplication - #282

I want to first acknowledge this with thanks -- it's quite impressive to be near GMP performance!

However to set expectations: this is also a large PR, and I will need some studying to understand what's going on, so it may take me a while to find time to review this.

Thanks. I'd like to note a few changes in the latest commit, which brings 10-15% performance gains and reduced memory footprint:

• The updated code uses decimation-in-frequency (DIF) in forward transform and decimation-in-time (DIT) in inverse transform.
• The prime numbers are replaced with bigger primes to enable tighter bit-packing in some cases that occur periodically.
  • P1 = 14_259_017_916_245_606_401, Max NTT length = 2^22 * 3^21 * 5^2 = 1_096_847_532_018_892_800
  • P2 = 17_984_575_660_032_000_001, Max NTT length = 2^19 * 3^17 * 5^6 = 1_057_916_215_296_000_000
  • P3 = 17_995_154_822_184_960_001, Max NTT length = 2^17 * 3^22 * 5^4 = 2_570_736_403_169_280_000
  • P2 and P3 are used for the two-prime NTT, whereas the three-prime NTT uses all three.

Does the chart shows that current algorithm is faster than GMP? That's impressive.

I ran benchmark fib_hex 100m from https://github.com/tczajka/bigint-benchmark-rs on this PR and it made num-bigint twice faster than malachite, slightly faster than gmp and 12x faster than itself.

Obviously the performance here is impressive, as I said before. However, as I finally started to try reviewing this, I have a few high level objections.

• The long commit history is not something I would want to merge. If there are discrete logical steps in implementing this feature, we can preserve some of that, but not a bunch of "update", "improve", and merges.
• The Montgomery implementation should share src/biguint/monty.rs, extending that if needed.
• So much unsafe code is uncharacteristic for this crate, and I would much rather stick to safe code as much as possible, even if that means sacrificing some performance. Case in point, miri finds UB here:

test test_mul_overflow_2 ... error: Undefined Behavior: attempting a write access using <5072044> at alloc1162269[0x0], but that tag does not exist in the borrow stack for this location
   --> src/biguint/ntt.rs:230:13
    |
230 |             *out.add(i) = Arith::<P>::mreduce(v);
    |             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ this error occurs as part of an access at alloc1162269[0x0..0x8]
    |
    = help: this indicates a potential bug in the program: it performed an invalid operation, but the Stacked Borrows rules it violated are still experimental
    = help: see https://github.com/rust-lang/unsafe-code-guidelines/blob/master/wip/stacked-borrows.md for further information
help: <5072044> was created by a SharedReadWrite retag at offsets [0x30..0x4830]
   --> src/biguint/ntt.rs:518:27
    |
518 |         conv_base::<P>(g, (&mut x[i..]).as_mut_ptr(), (&mut y[i..]).as_mut_ptr(), tf_current);
    |                           ^^^^^^^^^^^^^^^^^^^^^^^^^^
    = note: this is on thread `test_mul_overfl`
    = note: stack backtrace:
            0: num_bigint::biguint::ntt::conv_base::<17984575660032000001>
                at src/biguint/ntt.rs:230:13: 230:49
            1: num_bigint::biguint::ntt::conv::<17984575660032000001>
                at src/biguint/ntt.rs:518:9: 518:94
            2: num_bigint::biguint::ntt::mac3_two_primes
                at src/biguint/ntt.rs:604:5: 604:98
            3: num_bigint::biguint::ntt::mac3_u64
                at src/biguint/ntt.rs:760:9: 760:41
            4: num_bigint::biguint::ntt::mac3
                at src/biguint/ntt.rs:772:5: 772:24
            5: num_bigint::biguint::multiplication::mac3
                at src/biguint/multiplication.rs:358:9: 358:29
            6: num_bigint::biguint::multiplication::mul3
                at src/biguint/multiplication.rs:366:5: 366:31
            7: num_bigint::biguint::multiplication::<impl std::ops::Mul for &num_bigint::BigUint>::mul
                at src/biguint/multiplication.rs:428:31: 428:41
            8: test_mul_overflow_2
                at tests/biguint.rs:1075:20: 1075:27
            9: test_mul_overflow_2::{closure#0}
                at tests/biguint.rs:1068:25: 1068:25

(that particular one does pass under tree-borrows though)

Thanks for your review and comments! I’ll try to address the issues, including removing the unsafe code and cleaning up the commit history, even if that comes at the cost of a small performance regression.

Integrating the Montgomery reduction code with src/monty.rs is a bit more complicated, though:

• The NTT code currently uses a 64-bit modulus regardless of the target architecture, whereas src/monty.rs uses BigDigit. Introducing generics could potentially address this mismatch.
• The Montgomery operations in this PR were hand-optimized by inspecting the LLVM output on Godbolt and benchmarking promising variants. They are tailored specifically to the operations used by the NTT routine, including fused operations and the careful omission of overflow checks where they are mathematically unnecessary. As a result, they serve a somewhat different purpose from the code in src/monty.rs, and I’m not yet sure how they could be cleanly integrated there.

We can treat the monty-consolidation as "nice to have". However, any specific tailoring needs comments, lest that work be undone by a later contributor or maintainer. I'm skeptical that those optimizations wouldn't be useful to the other monty use as well though...

I’ve made the suggested improvements:

• Cleaned up the commit history.
  The original branch remains available here for reference.
• Rebased the branch onto the latest main.
• Replaced the unsafe routines with safe equivalents. Benchmarks showed no performance regression.
• Added comments explaining why the NTT routine requires its own Montgomery arithmetic routines.

I also made a small improvement to the planner based on exhaustive Criterion.rs benchmarks covering operand sizes of up to 300 million bits each. This reduced execution time by 1.2% on average.

Some of this work was completed with assistance from ChatGPT Codex.

Thanks again for the detailed review. Please let me know if there are any remaining issues.

Hi, thanks for the changes! I’ve been reviewing this and verified via benchmarks that using fixed-length slices does not affect performance; I also checked that the generated assembly is unchanged. I’ve also been examining the planner and gathering additional data.

For prime selection, we should balance two considerations: the modulus should be large to increase the packing density, while also having as many factors of 2, 3, and 5 as possible. Currently, the primes are selected by exhaustive search so that they have sufficient factors of 2, 3, and 5 while remaining close to 2^64.

The m7 check was for the base case, which is handled by naive quadratic multiplication, so we can include at most one factor of 7 in the length. This reduces the amount of padding required for certain input lengths.

Besides this, I’ve carefully checked the planner scheme and found that it can be somewhat simplified while retaining the same planning accuracy. I’ve been preparing those changes, and I should have communicated this earlier—sorry about that.

The updated planner is available here. Probably I could open a new PR.